#### Independent and combined associations of urinary arsenic exposure and serum sex steroid hormones among 6-19-year old children and adolescents in NHANES 2013–2016 ####
# 2023-3-30 14:16:05 进行复现论文paper数据
# 尿砷暴露的独立和联合相关性 6-19岁儿童和青少年血清性类固醇激素  选取的数据是ehanes 2013-2014 2015-2016年数据
# 文章摘要说了结论 关于 暴露于砷，无论是单独或混合物，可能改变性类固醇激素水平儿童和青少年。
# 这个是关于 gWQS 这个不熟悉啊,完蛋看都看不懂...
# https://blog.csdn.net/weixin_42812146/article/details/126192945 

# 2023年7月12日10:46:56 开始搞
# 里面有bayes 前面有数据方法支持
# 在周期里面 排除了 BMI, BMI分类, serum cotinine ,PIR 丢失


library(tidyverse)
library(gWQS)
library(ggplot2)
library(haven)
library(gtsummary)
library(bkmr)

  setwd('G:/BaiduNetdiskDownload/nhanes')

setwd("${path}")

# 数据写入本地进行减少内存使用
# write.csv(output619Arsenouslbx,file = 'F:/Rproject/r-language/paper2/output619Arsenouslbx.csv')
 output619Arsenouslbx<- read.csv(file = 'F:/Rproject/r-language/paper2/output619Arsenouslbx.csv')
# agerange317dropna <- read.csv('${path}source.csv')
output619Arsenouslbx<- read.csv(file = '${path}output619Arsenouslbx.csv')
#1069  人 文章是1063 暂且不管了 2023年7月14日14:01:55 


# 儿童:6 ~ 11岁(children)，青少年:12 ~ 19岁 (adolescents)

table1data <- output619Arsenouslbx

# tab1 数据 age race BMI ,六个月时间段(无法找到对应),抽血时间(无法找到对应), 贫困指数PIR ,可替宁暴露状态, 尿肌酐,血清性类固醇激素(TT ,E2,SHBG,FAI,TT/E2), 尿砷种类代谢物(总砷,亚砷酸,DMA,MMA)


table1data$Sex <- ifelse(table1data$RIAGENDR==1,'Male','Female')

table1data$race <- recode_factor(table1data$RIDRETH1, 
                            `1` = 'Mexican American',
                            `2` = 'Other Hispanic',
                            `3` = 'Non-Hispanic White',
                            `4` = 'Non-Hispanic Black',
                            `5` = 'Other race - including multi-racial'
                            
)
table1data$Cotinineexposurestatus<- ifelse(table1data$LBXCOT >0.015,'Exposed','Unexposed')

table1data$Age <- table1data$RIDAGEYR
table1data$BMI <- table1data$BMXBMI
# 尿肌酐
table1data$Urinarycreatinine <- table1data$URXUCR
# Total testosterone (TT ) 总睾酮 TST	lbxtst,  estradiol (E2) 雌二醇 Laboratory	TST	lbxest , sex hormone-binding globulin   性激素结合球蛋白(SHBG) Laboratory	TST	lbxshbg
# 检测了三种尿砷代谢物  亚砷酸 (Arsenous acid UAS	urxuas3)、DMA (二甲基胂酸 UAS urxudma)、MMA(一甲基胂酸 UAS	urxumma)
# TT
table1data$TT <- table1data$LBXTST
#E2
table1data$E2 <- table1data$LBXEST
#SHBG
table1data$SHBG <- table1data$LBXSHBG
#FAI  the free androgen index (FAI) generated by TT/SHBG
table1data$FAI <- table1data$TT/table1data$SHBG
# TT/E2
table1data$TTE2 <- table1data$TT/table1data$E2

# urxuas Arsenic, total

table1data$Arsenictotal <- table1data$URXUAS

#Arsenous acid
table1data$Arsenousacid <- table1data$URXUAS3
#dma
table1data$DMA <- table1data$URXUDMA
#mma
table1data$MMA <- table1data$URXUMMA
table1data$PIR <- table1data$INDFMPIR

table1data$childgroup <- ifelse(table1data$RIDAGEYR<=11,'Children','Adolescents')
colnames(table1data)

cols <-colnames(table1data)
coluse <-cols[c(100:116)]
coluse
table1data<-table1data[,coluse]


table1datamale <- table1data[which(table1data$Sex=='Male'),]
table1dataFemale <- table1data[which(table1data$Sex=='Female'),]

table1datademoAll <- table1data
write.csv(table1datademoAll,'F:/Rproject/r-language/paper2/output619ArsenouslbxUse.csv')

# 2023年7月18日08:46:58 暂且没搜到多分组 直接两个分组进行合并吧
# tbl_summary(table1datademo,label = list(Age ~ "Age (year)"))
colnames(table1datademoAll)
tab1All<-tbl_summary(table1datademoAll,label = list(Age ~ "Age (year)",
                                                    race~"Race",
                                                    Cotinineexposurestatus~"Cotinine exposure status",
                                                    BMI~"BMI (kg/m2)",
                                                    Urinarycreatinine~"Urinary creatinine (mg/dL)",
                                                    TT~"TT (ng/dL))",
                                                    E2~"E2 (pg/mL)",
                                                    SHBG~"SHBG (nmol/L)",
                                                    TT~"TT (ng/dL))",
                                                    TTE2~"TT/E2",
                                                    TT~"TT (ng/dL))",
                                                    Arsenictotal~"Arsenic, total",
                                                    Arsenousacid~"Arsenous acid" ))
tab1male<-tbl_summary(table1datamale,by = childgroup,label = list(Age ~ "Age (year)",
                                                                      race~"Race",
                                                                      Cotinineexposurestatus~"Cotinine exposure status",
                                                                      BMI~"BMI (kg/m2)",
                                                                      Urinarycreatinine~"Urinary creatinine (mg/dL)",
                                                                      TT~"TT (ng/dL))",
                                                                      E2~"E2 (pg/mL)",
                                                                      SHBG~"SHBG (nmol/L)",
                                                                      TT~"TT (ng/dL))",
                                                                      TTE2~"TT/E2",
                                                                      TT~"TT (ng/dL))",
                                                                      Arsenictotal~"Arsenic, total",
                                                                      Arsenousacid~"Arsenous acid" )
                      ) %>%modify_spanning_header(all_stat_cols() ~ "**Male**")
tab1female<-tbl_summary(table1dataFemale,by = childgroup,label = list(Age ~ "Age (year)",
                                                                          race~"Race",
                                                                          Cotinineexposurestatus~"Cotinine exposure status",
                                                                          BMI~"BMI (kg/m2)",
                                                                          Urinarycreatinine~"Urinary creatinine (mg/dL)",
                                                                          TT~"TT (ng/dL))",
                                                                          E2~"E2 (pg/mL)",
                                                                          SHBG~"SHBG (nmol/L)",
                                                                          TT~"TT (ng/dL))",
                                                                          TTE2~"TT/E2",
                                                                          TT~"TT (ng/dL))",
                                                                          Arsenictotal~"Arsenic, total",
                                                                          Arsenousacid~"Arsenous acid" )) %>%modify_spanning_header(all_stat_cols() ~ "**FeMale**")

tab1all<-tbl_merge(tbls = list(tab1All,tab1female,tab1male),tab_spanner =c('ALL','Female','Male') )

tab1all%>%as_flex_table()%>%flextable::save_as_html(path = 'paper2tab1.html')

#tab2 多种砷代谢物 (亚砷酸、DMA、MMA) 对性类固醇激素水平影响((TT、E2、SHBG、FAI、TT/E2)
# Total testosterone (TT ) 总睾酮 TST	lbxtst,  estradiol (E2) 雌二醇 Laboratory	TST	lbxest , sex hormone-binding globulin   性激素结合球蛋白(SHBG) Laboratory	TST	lbxshbg
# 检测了三种尿砷代谢物  亚砷酸 (Arsenous acid UAS	urxuas3)、DMA (二甲基胂酸 UAS urxudma)、MMA(一甲基胂酸 UAS	urxumma)
colnames(table1datademoAll)

# 首先对数值进行ln转换 所有砷种类和性激素变量均进行了ln转换

table1datademoAllParse <- table1datademoAll
#进行激素类转换
table1datademoAllParse$TT<- log(table1datademoAllParse$TT)
table1datademoAllParse$E2<- log(table1datademoAllParse$E2)
table1datademoAllParse$SHBG<- log(table1datademoAllParse$SHBG)
table1datademoAllParse$FAI<- log(table1datademoAllParse$FAI)
#进行代谢物转换
table1datademoAllParse$Arsenictotal<- log(table1datademoAllParse$Arsenictotal)
table1datademoAllParse$Arsenousacid<- log(table1datademoAllParse$Arsenousacid)
table1datademoAllParse$DMA<- log(table1datademoAllParse$DMA)
table1datademoAllParse$MMA<- log(table1datademoAllParse$MMA)

# 因为无法找到其他两个协变量 用已知的协变量处理

# glm
#at
ArsenictotalGLM <- glm(Arsenictotal ~ Age+Sex+race+BMI+PIR+Cotinineexposurestatus+TT+E2+SHBG+FAI+TTE2,
                       data = table1datademoAllParse)



tab2regat<-tbl_regression(ArsenictotalGLM,include =c(TT,E2,SHBG,FAI,TTE2),exponentiate = T)

#acid
ArsenousacidGLM <- glm(Arsenousacid ~ Age+Sex+race+BMI+PIR+Cotinineexposurestatus+TT+E2+SHBG+FAI+TTE2,
                       data = table1datademoAllParse)



tab2regacid<-tbl_regression(ArsenousacidGLM,include =c(TT,E2,SHBG,FAI,TTE2),exponentiate = T)
#DMA

dmaGLM <- glm(DMA ~ Age+Sex+race+BMI+PIR+Cotinineexposurestatus+TT+E2+SHBG+FAI+TTE2,
                       data = table1datademoAllParse)



tab2regDMA<-tbl_regression(dmaGLM,include =c(TT,E2,SHBG,FAI,TTE2),exponentiate = T)

# MMA
MMAGLM <- glm(MMA ~ Age+Sex+race+BMI+PIR+Cotinineexposurestatus+TT+E2+SHBG+FAI+TTE2,
                       data = table1datademoAllParse)



tab2regMMA<-tbl_regression(MMAGLM,include =c(TT,E2,SHBG,FAI,TTE2),exponentiate = T)

paper2tab2<-tbl_stack(tbls = list(tab2regat,tab2regacid,tab2regDMA,tab2regMMA),group_header = c('ln-Arsenic, total','ln-Arsenous Acid','ln-DMA','ln-MMA'))

paper2tab2%>%as_flex_table()%>%flextable::save_as_html(path = 'paper2tab2.html')
# 结束合并




# 去了解wqs


toxic_chems <- names(table1datademoAllParse)[13:15]
names(table1datademoAllParse)

# we run the model and save the results in the variable "results"
# monomethylarsonic acid. (A), (B), (C) and (D) were the weight
# of each arsenic metabolite in the positive WQS model. 
#(A) TT: total testosterone;
#(B) E2: estradiol; 
#(C) FAI: free androgen index, was calculated as [TT (ng/dL)] / [SHBG(nmol/L)];
#(D) TT/E2 was calculated as [TT (ng/dL)] / [E2 (pg/mL)]
#A
resultsTTposA <- gwqs(TT ~ wqs, # 因子分析数据
                mix_name = toxic_chems, # 混合因素
                data = table1datademoAllParse, #数据来源
                q = 100, 
                validation = 0.6,# 训练验证区分大小
                b = 200, # bootstrap样本数量
                b1_pos = TRUE,  # 假设 β1为正（b1_pos = TRUE）
                b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                family = "gaussian",
                seed = 1000,  # 设置种子数 用于可以获得重复结果
                plots = TRUE,
                tables = TRUE
                )

Cairo::CairoTIFF(file="TTtotaltestosterone.tiff", width=8, height=8,units="in",dpi=150)
TTtotaltestosterone<- gwqs_barplot(resultsTTposA)
TTtotaltestosterone


#(A) TT: total testosterone;
#(B) E2: estradiol; 
#(C) FAI: free androgen index, was calculated as [TT (ng/dL)] / [SHBG(nmol/L)];
#(D) TT/E2 was calculated as [TT (ng/dL)] / [E2 (pg/mL)]
#B
resultsE2posB <- gwqs(E2 ~ wqs, # 因子分析数据
                     mix_name = toxic_chems, # 混合因素
                     data = table1datademoAllParse, #数据来源
                     q = 100, 
                     validation = 0.6,# 训练验证区分大小
                     b = 200, # bootstrap样本数量
                     b1_pos = TRUE,  # 假设 β1为正（b1_pos = TRUE）
                     b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                     family = "gaussian",
                     seed = 1000,  # 设置种子数 用于可以获得重复结果
                     plots = TRUE,
                     tables = TRUE
)

Cairo::CairoTIFF(file="E2.tiff", width=8, height=8,units="in",dpi=150)
E2<-gwqs_barplot(resultsE2posB)
E2

#C
resultsFAIposC <- gwqs(FAI ~ wqs, # 因子分析数据
                      mix_name = toxic_chems, # 混合因素
                      data = table1datademoAllParse, #数据来源
                      q = 100, 
                      validation = 0.6,# 训练验证区分大小
                      b = 200, # bootstrap样本数量
                      b1_pos = TRUE,  # 假设 β1为正（b1_pos = TRUE）
                      b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                      family = "gaussian",
                      seed = 1000,  # 设置种子数 用于可以获得重复结果
                      plots = TRUE,
                      tables = TRUE
)

Cairo::CairoTIFF(file="FAI.tiff", width=8, height=8,units="in",dpi=150)
FAI<-gwqs_barplot(resultsFAIposC)
FAI

#D
resultsTTE2posd <- gwqs(TTE2 ~ wqs, # 因子分析数据
                       mix_name = toxic_chems, # 混合因素
                       data = table1datademoAllParse, #数据来源
                       q = 100, 
                       validation = 0.6,# 训练验证区分大小
                       b = 200, # bootstrap样本数量
                       b1_pos = TRUE,  # 假设 β1为正（b1_pos = TRUE）
                       b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                       family = "gaussian",
                       seed = 1000,  # 设置种子数 用于可以获得重复结果
                       plots = TRUE,
                       tables = TRUE
)

Cairo::CairoTIFF(file="TTE2.tiff", width=8, height=8,units="in",dpi=150)
TTE2<-gwqs_barplot(resultsTTE2posd)
TTE2

#  (E), (F), (G), (H), and (I) were the weight of each arsenic metabolite in the negative WQS model.
#(E) TT: total testosterone; 
#(F) E2: estradiol;
#(G) SHBG: sex hormone binding globulin;
#(H) FAI: free androgen index, was calculated as [TT (ng/dL)] / [SHBG (nmol/L)];
#(I) TT/E2 was calculated as [TT (ng/dL)] / [E2 (pg/mL)].

#E
resultsTTnegE <- gwqs(TT ~ wqs, # 因子分析数据
                        mix_name = toxic_chems, # 混合因素
                        data = table1datademoAllParse, #数据来源
                        q = 100, 
                        validation = 0.6,# 训练验证区分大小
                        b = 200, # bootstrap样本数量
                        b1_pos = FALSE,  # 假设 β1为正（b1_pos = TRUE）
                        b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                        family = "gaussian",
                        seed = 1000,  # 设置种子数 用于可以获得重复结果
                        plots = TRUE,
                        tables = TRUE
)

Cairo::CairoTIFF(file="resultsTTnegE.tiff", width=8, height=8,units="in",dpi=150)
postTTnegE<-gwqs_barplot(resultsTTnegE)
postTTnegE

#F
resultsE2negF <- gwqs(E2 ~ wqs, # 因子分析数据
                        mix_name = toxic_chems, # 混合因素
                        data = table1datademoAllParse, #数据来源
                        q = 100, 
                        validation = 0.6,# 训练验证区分大小
                        b = 200, # bootstrap样本数量
                        b1_pos = FALSE,  # 假设 β1为正（b1_pos = TRUE）
                        b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                        family = "gaussian",
                        seed = 1000,  # 设置种子数 用于可以获得重复结果
                        plots = TRUE,
                        tables = TRUE
)

Cairo::CairoTIFF(file="resultsE2negF.tiff", width=8, height=8,units="in",dpi=150)

E2negF<-gwqs_barplot(resultsE2negF)
E2negF
#(E) TT: total testosterone; 
#(F) E2: estradiol;
#(G) SHBG: sex hormone binding globulin;
#(H) FAI: free androgen index, was calculated as [TT (ng/dL)] / [SHBG (nmol/L)];
#(I) TT/E2 was calculated as [TT (ng/dL)] / [E2 (pg/mL)].
#G
resultsSHBGnegG <- gwqs(SHBG ~ wqs, # 因子分析数据
                        mix_name = toxic_chems, # 混合因素
                        data = table1datademoAllParse, #数据来源
                        q = 100, 
                        validation = 0.6,# 训练验证区分大小
                        b = 200, # bootstrap样本数量
                        b1_pos = FALSE,  # 假设 β1为正（b1_pos = TRUE）
                        b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                        family = "gaussian",
                        seed = 1000,  # 设置种子数 用于可以获得重复结果
                        plots = TRUE,
                        tables = TRUE
)

Cairo::CairoTIFF(file="resultsSHBGnegG.tiff", width=8, height=8,units="in",dpi=150)
SHBGnegG<-gwqs_barplot(resultsSHBGnegG)
SHBGnegG
#H
resultsFAInegH <- gwqs(FAI ~ wqs, # 因子分析数据
                        mix_name = toxic_chems, # 混合因素
                        data = table1datademoAllParse, #数据来源
                        q = 100, 
                        validation = 0.6,# 训练验证区分大小
                        b = 200, # bootstrap样本数量
                        b1_pos = FALSE,  # 假设 β1为正（b1_pos = TRUE）
                        b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                        family = "gaussian",
                        seed = 1000,  # 设置种子数 用于可以获得重复结果
                        plots = TRUE,
                        tables = TRUE
)

Cairo::CairoTIFF(file="resultsFAInegH.tiff", width=8, height=8,units="in",dpi=150)
FAInegH<-gwqs_barplot(resultsFAInegH)
FAInegH
#I
resultsTTE2negi <- gwqs(TTE2 ~ wqs, # 因子分析数据
                        mix_name = toxic_chems, # 混合因素
                        data = table1datademoAllParse, #数据来源
                        q = 100, 
                        validation = 0.6,# 训练验证区分大小
                        b = 200, # bootstrap样本数量
                        b1_pos = FALSE,  # 假设 β1为正（b1_pos = TRUE）
                        b1_constr = TRUE,#  设置该参数为假（b1_pos = FALSE）来测试负相关
                        family = "gaussian",
                        seed = 1000,  # 设置种子数 用于可以获得重复结果
                        plots = TRUE,
                        tables = TRUE
)
#存储图片
Cairo::CairoTIFF(file="resultsTTE2negd.tiff", width=8, height=8,units="in",dpi=150)
TTE2negi<-gwqs_barplot(resultsTTE2negi)
TTE2negi

# 测试贝叶斯 A tt
# 单个砷代谢物与每种性类固醇激素的暴露-反应关系函数
# (A) TT:总睾酮;
#(B)E2:雌二醇;
#(C) SHBG:性激素结合球蛋白;
#(D) FAI:游离雄激素指数，计算公式为[TT (ng/dL)] / [SHBG (nmol/L)];
#(E) TT/E2计算为[TT (ng/dL)] / [E2(pg/mL)]
demoS4 <-table1datademoAllParse
demoS4<- demoS4[which(!is.na(demoS4$DMA)),]
# 数据200 1000的话太慢了
demoS4 <- demoS4[c(1:200),]
colnames(demoS4)
Z <- demoS4[,c('Arsenousacid','DMA','MMA')]

#A
set.seed(1000)
yA <- demoS4$TT
fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
risks.overallA = OverallRiskSummaries(fit=fitkmA,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="TTbayes.tiff", width=8, height=8,units="in",dpi=150)
TTbayes<-ggplot(risks.overallA,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+ggtitle("TT")+
  geom_pointrange()



TTbayes

# B
yB <- demoS4$E2
set.seed(1000)
fitkmB <- kmbayes(y =yB, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
risks.overallB = OverallRiskSummaries(fit=fitkmB,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="E2bayes.tiff", width=8, height=8,units="in",dpi=150)
E2bayes<-ggplot(risks.overallB,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+ggtitle("E2")+
  geom_pointrange()


E2bayes

# c
yC <- demoS4$SHBG

set.seed(1000)
fitkmC <- kmbayes(y =yC, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
risks.overallC = OverallRiskSummaries(fit=fitkmC,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="SHBGbayes.tiff", width=8, height=8,units="in",dpi=150)
SHBGbayes<-ggplot(risks.overallC,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+ggtitle("SHBG")+
  geom_pointrange()


SHBGbayes

# D
yD <- demoS4$FAI
set.seed(1000)
fitkmD <- kmbayes(y =yD, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
risks.overallD = OverallRiskSummaries(fit=fitkmD,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="FAIbayes.tiff", width=8, height=8,units="in",dpi=150)
FAIbayes<-ggplot(risks.overallD,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+ggtitle("FAI")+
  geom_pointrange()

FAIbayes

# E
yE <- demoS4$TTE2
set.seed(1000)
fitkmE <- kmbayes(y =yE, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
risks.overallE = OverallRiskSummaries(fit=fitkmE,qs=seq(0.25,0.75,by=0.05),q.fixed = 0.5)
Cairo::CairoTIFF(file="TTE2bayes.tiff", width=8, height=8,units="in",dpi=150)
TTE2bayes<-ggplot(risks.overallE,aes(quantile,est,ymin=est-1.96*sd,ymax=est+1.96*sd))+
  geom_hline(yintercept = 0,lty=2,col='brown')+ggtitle("TT/E2")+
  geom_pointrange()
TTE2bayes

##########################################################结束

# 混合暴露与性类固醇激素的联合关系
#A
pred.resp.univarA <- PredictorResponseUnivar(fit = fitkmA)
Cairo::CairoTIFF(file="A1bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univarA, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)


#B
pred.resp.univarB <- PredictorResponseUnivar(fit = fitkmB)
Cairo::CairoTIFF(file="B1bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univarB, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)

#C
pred.resp.univarC <- PredictorResponseUnivar(fit = fitkmC)
Cairo::CairoTIFF(file="C1bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univarC, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)

#D
pred.resp.univarD <- PredictorResponseUnivar(fit = fitkmD)
Cairo::CairoTIFF(file="D1bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univarD, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)

#E
pred.resp.univarE <- PredictorResponseUnivar(fit = fitkmE)
Cairo::CairoTIFF(file="E1bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.univarA, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(expos)") +
  xlab("Phthalate metabolites") +
  xlim(-2,6)


#砷代谢物和性类固醇激素的双变量暴露-反应函数
# 最后图

set.seed(1000)
#A
# fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
pred.resp.bivarA <- PredictorResponseBivar(fit =fitkmA , min.plot.dist = 1)
pred.resp.bivar.levelsA <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivarA, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
Cairo::CairoTIFF(file="A2bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levelsA, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")


#B
# fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
pred.resp.bivarB <- PredictorResponseBivar(fit =fitkmB , min.plot.dist = 1)
pred.resp.bivar.levelsB <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivarB, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
Cairo::CairoTIFF(file="B2bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levelsB, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
#C
# fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
pred.resp.bivarC <- PredictorResponseBivar(fit =fitkmC , min.plot.dist = 1)
pred.resp.bivar.levelsC <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivarC, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
Cairo::CairoTIFF(file="c2bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levelsC, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
#D
# fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
pred.resp.bivarD <- PredictorResponseBivar(fit =fitkmD , min.plot.dist = 1)
pred.resp.bivar.levelsD <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivarD, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
Cairo::CairoTIFF(file="D2bayes.tiff", width=8, height=8,units="in",dpi=150)
ggplot(pred.resp.bivar.levelsD, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
#E
# fitkmA <- kmbayes(y =yA, Z = Z, X = NULL, iter = 1000, verbose = FALSE, varsel = TRUE,family = 'gaussian',est.h = TRUE)
pred.resp.bivarE <- PredictorResponseBivar(fit =fitkmE , min.plot.dist = 1)
pred.resp.bivar.levelsE <- PredictorResponseBivarLevels(  
  pred.resp.df = pred.resp.bivarE, 
  Z = Z, qs = c(0.1, 0.5, 0.9))
Cairo::CairoTIFF(file="EE2bayes.tiff", width=8, height=8,units="in",dpi=150)
EE2bayes<-ggplot(pred.resp.bivar.levelsE, aes(z1, est)) + 
  geom_smooth(aes(col = quantile), stat = "identity") + 
  facet_grid(variable2 ~ variable1) +
  ggtitle("h(expos1 | quantiles of expos2)") +
  xlab("expos1")
EE2bayes



